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SUMMARY 

A finite difi'erence upwind discretization scheme in two dimensions is presented in detail for 
the transient simulation of the highly coupled non-linear partial differential equations of the full 
hydrodynamic model, providing thereby a practical engineering tool for improved charge carrier 
transport simulations at high electric fields and frequencies. The discretization scheme preserves 
the conservation and transportive properties of the equations. The hydrodynamic model is able to 
describe inertia effects which play an increasing role in different fields of micro- and optoelectronics, 
where simplified charge transport models like the drift-diffusion model and the energy balance model 
are no longer applicable. Results of extensive numerical simulations are shown for a two-dimensional 
MESFET device. A comparison of the hydrodynamic model to the commonly used energy balance 
model is given and the accuracy of the results is discussed. 

KEY WORDS: Semiconductor device modeling; charge transport models; hydrodynamic model; 
upwind discretization; submicron devices; hot electrons; velocity overshoot, Monte Carlo methods, 
MESFETs. 



1. INTRODUCTION 

There is a growing interest in extended charge transport models for semiconductor devices. 
Our paper emerges from the the fact that today's submicron semiconductor devices like e.g. 
MESFETs and HEMTs are operated under strong electric fields and at high frequencies. 
Information transmission using an electromagnetic wave at very high frequencies will have a 
direct impact on how we design active and passive components in different fields of micro- 
and optoelectronics. In such cases, quasi-static semiconductor device models like the energy 
balance model (EBM) are no longer adequate. Especially in GaAs and related materials used 
for high-speed device design, inertia effects play an important role since the momentum and 
energy relaxation times of the electron gas are close to the picosecond range. 
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The most elaborate and practicable approach for the description of charge transport in 
semiconductors used for device simulation would be the Monte Carlo (MC) method llj. 
The advantage of this technique is a complete picture of carrier dynamics with reference 
to microscopic material parameters, e.g. effective masses and scattering parameters. But the 
method must be still considered as very time consuming and hence too uneconomical to be 
used by device designers. 

Besides the simplest concept which is the traditional drift-diffusion model (DDM), there 
is a much more rigorous approach to the problem, namely the so-called hydrodynamic model 
(HDM). A simplified version of the HDM is the EBM, which is implemented in most of today's 
commercial semiconductor device simulators. The HDM makes use of electron temperature for 
the description of charge transport and takes inertia effects into account as well. Starting 
from the Boltzmann equation, Blotekjaer j2| and many others presented a derivation of such 
equations under the assumption of parabolic band structures. Especially for silicon, satisfactory 
results were obtained this way 131. But the results have often been unsatisfactory when MC 
models based on a nonparabolic structure were compared to HDM results based on an empirical 
choice of model parameters. Therefore, it is quite natural to improve the HDM by incorporating 
energy-dependent relaxation times and effective masses obtained from MC bulk simulations 

^ . This is a strategy we pursue in this paper. 

In the first part of this work, we give a short definition of the hydrodynamic model for GaAs. 
We emphasize that a thorough analysis of the physical features of the charge carrier transport 
models is the basis for a clear understanding of their limits of applicability. Then we give a 
simple discretization scheme for the full hydrodynamic model in two dimensions, which will be 
applied to a GaAs MESFET structure in the last part. There, we compare the HDM results 
also to results obtained from the simpler EBM, and we investigate the influence of the grid 
resolution on the results. 



2. THE HYDRODYNAMIC MODEL FOR GaAs FETs 
2.1. Definition of the model 

Our active device model is based on the single-gas hydrodynamic equations. This is a 
simplification of the two-valley hydrodynamic equations, since strictly speaking, in GaAs and 
other semiconductors with similar band structure like InP, there exists an electron gas with 
different thermal distribution function for each conduction band valley (i.e. the F- and L- 
valleys). The equations for each valley are, however, coupled through collision terms since 
electrons can scatter between two different valleys. The corresponding relaxation rates may 
be of the order of a picosecond and are therefore relatively large. This is the main drawback 
of the single-gas approximation, and it would be desirable to implement at least a two- valley 
hydrodynamic model. Reliable extensive two-valley simulations have been performed only for 
the one-dimensional case so far due to the large amount of equations and parameters involved 
in such a model. A hydrodynamic two-valley simulation of GaAs MESFETs is the subject of 
one of our forthcoming papers. The HDM equations consist of the continuity equation 

^ + V(™) = (1) 
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for negligible charge carrier generation and recombination, the momentum balance equation 
given by 

^+(yp)v + {py)v = -enE-y{nkT)-^ (2) 

Ot Tr, 



or alternatively (only for the x-component) 

d{m*ijjo)nVx) 



d{nkT) m*{u))nvx 
-qnEx ^- , (3) 



and the energy balance equation 



dt 

dx '^p(^) 
dco 



-envE - V{nkTv) - V(-kVT) ^— , (4) 

where n, u [u = to/n), and v are the electron density, the electron energy density (average 
electron energy) and the electron drift velocity, respectively, is the x-component of the 
electron drift velocity and p = m*nv the momentum density. Corresponding equations are 
valid for the y- (and z-) components. T is the electron temperature, e > the elemental 
charge and k Boltzmann's constant, ujq = f fc^L is the average thermal equilibrium energy of 
electrons, where Tl is the lattice temperature. The EBM uses only a simplified energy balance 
equation (see below). The electronic current density J inside the active device is given by 

J = —env , (5) 

the total current density is 

? ^ dE 

Jtot = -env + eo€r^ ■ (6) 
ot 

The momentum relaxation time Tp{ljj) is related to the mobility of the electrons via /x(ZZJ) = 
{e/m*{cJ))Tp{u}), and the energy relaxation time ra,(w) describes the exchange of energy 
between the heated electron gas and the lattice. Tj, and and the effective electron mass m* 
are assumed to be functions of the mean electron energy. We performed steady-state Monte 
Carlo simulations in order to get the correct values for these parameters. 

The hydrodynamic equations, together with Poisson's equation {N^ (— N^) is the number 
of (ionized) donors) 

A(/)=-V£ = —{N+-n) (7) 

eoCr 

form a complete set of equations that can be used to solve for the electron density, velocity, 
energy and electric field for given boundary conditions, if we use a closing relation for the 
mean electron energy ZZJ, the electron temperature T and velocity v: 

cZJ = ^m**{uy + ^kT + (3l{'^)AEtl ■ (8) 
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The reason for the double index of the electron mass will soon become clear. The last term in 
eq. ((SJ accounts for the fact that a minimum energy AE'fl = 0.29 eV is necessary to excite an 
electron to an upper conduction band. (3^ is the relative fraction of electrons in the L-band for 
the stationary homogeneous case. The term /3i(a;)A£'rL is often neglected in the literature, 
but this may lead to an overestimation of the electron temperature of more than 1000 K at 
high energies. 

2.2. Remarks on the single gas approximation 

We point out again the important fact that we are using a single-gas approximation for the 
hydrodynamic model. This means that the closing relation eq. ijSJl is a crude approximation 
which allows the calculation of the electron temperature from the total electron energy and 
electron drift velocity. Some authors neglect also the influence of the electron velocity on the 
temperature [3 EI by directly relating the electron temperature T to the average electron 
energy uJ from stationary MC simulations. 

The transition from the two-gas model to the single-gas approximation has to be done 
carefully; therefore, we present here a short discussion of the problem which is usually not 
mentioned in the literature. We assume parabolic T— and L— valleys for the sake of brevity, 
but in our simulations, we took the effect of the non-parabolicity of the energy bands into 
account by using the Kane model [7] , which generalizes the parabolic relation for the electron 
energy and electron crystal momentum hk 

Ek = = 7(fc) (9) 

to 

Ekil - aEk) ^ -/{k) , (10) 

with the non-parabolicity coefficient a which has different values for the different energy bands, 
mgfj is the effective electron mass at the bottom of the energy band under consideration. Very 
often, an energy-dependent electron mass is defined via 



J_ ^ 1 d^Ek 

1 1 



mgff (1 + 4a7(fc))3/2 



(11) 
(12) 



i.e. the electron mass increases with growing energy. The electron mass ttle must not be 
confused with the energy-dependent electron mass used in the hydrodynamic model. In that 
case, the masses are a kind of average masses depending on the average electron energy. 

There is even a further aspect related to the notion of the electron mass. The crystal velocity 
of an electron is given for spherical bands by 

IdEu 



h dk ■ 

A short calculation for the Kane model shows that 

hk 1 



(13) 



(14) 
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which implies that crystal velocity v and crystal momentum p = hk are related by 



p = nipV , rrip = + 4:a^{k)m^ff , (15) 

i.e. a different definition of the energy-dependent electron mass applies if the electron velocity 
is calculated from the crystal momentum. 

In the single particle two-band MC simulations, the random walk of an electron inside the 
semiconductor material is monitored over a sufficiently long time. As a result, the probability 
Pr that an electron resides in the F— valley is obtained as a function of the applied constant 
homogeneous electric field i? or as a function of the mean electron energy oJ, and the probability 
of finding the electron in an upper L— valley is then Pl — 1 ~ /3r- Also the values for for the 
average electron velocities vr and in the different valleys are obtained as well. Then it is 
reasonable to define the average electron velocity by 

V = Prvr + Plvl ■ (16) 

The average electron momentum p is given by 

p = m*v = m-cPrVv + rriLPLVL , (17) 

hence the (energy-dependent) electron mass which must be used in the hydrodynamic model 
in order to relate average electron velocity and electron momentum is calculated from 

But we suggest that a different mass m**(ZU) should be used for the calculation of the average 
kinetic electron energy in eq. ©. We identify 

lm**v^ = i/3rmrwr + \pLmLvl , (19) 



and therefore 

{Pvvt+PlvlY 



(20) 



m* and m** are not distinguished in the literature. It is tempting to use a naive definition for 
the electron mass 

m = Prmr + PLmL ■ (21) 

It is interesting that data in the literature for the energy-dependent electron mass are usually 
in better agreement with this definition. It is clear that the HDM is still an approximative 
description of charge carrier dynamics inside a semiconductor, and the different assumptions 
which are inherent in the derivation of the model may already cause larger errors in the 
simulation results than using only one mass. Therefore we do not claim that our discussion 
leads to improved simulation results, but it rather shows the limits of the frequently used 
model and points out that it is mandatory to maintain always a highest possible degree of 
consistency. 

Figure 1 shows m*, to** and to normalized to the free electron mass toq as functions of the 
average electron energy for a GaAs lattice temperature — 300 K and a low doping density 
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Figure 1. Energy-dependent electron masses used in the hydrodynamic simulations. 



Nd — lO^^'cm"'^. The results from MC simulations were smoothed by a polynomial fit and 
transferred afterwards into the hydrodynamic simulation program. We used mr = O.OGTtoq 
and rriL = O.SStoq. For high energies when nearly all electrons are in the upper bands, the 
electron mass even exceeds the value uil due to the non-parabolicity of the energy bands. 

Figure 2 shows the average energy ZJ of an electron in a constant homogeneous electric field 
E for GaAs. For each data point, the electron was scattered one million times (including so- 
called self-scattering), therefore the resulting curve is already quite smooth. It is clear that 
MC simulations deliver no data for average electron energies below ZJ < lUq = 38 meV, since 
the mean electron energy ZJ has this value if no electric field is applied to the crystal, and ZJ 
grows for increasing electric field. This is no major problem, since the low energy region is of 
minor importance for the hydrodynamic simulation and the necessary data can be obtained 
from theoretical considerations 0]. 

In order to complete the set of data which is necessary for hydrodynamic simulations, the 
electron velocity and energy relaxation times are depicted in Figs. 3 and 4 for doping densities 
Nd = lO^^'cm"'^ and 2 • lO^^cm"^. The characteristic shape of the velocity curve can be 
explained by the fact that at high energies the electrons jump into the L— bands where the 
electrons have a lower mobility than in the F— band. 

Finally, we need an expression for the thermal conductivity of the electron gas, which is 
given by theoretical considerations 

^ = (5/2 + r)n^^^T . (22) 

e 

Several different choices for r can be found in the literature, and many authors 01 [S] even 
neglect heat conduction in their models. As a matter of fact, heat conduction does not 
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Figure 3. v-E curve for GaAs. 
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Figure 4. Energy relaxation times for GaAs (from a MC simulation, as in Figs. 1-3) 

influence the simulation results very much if r remains within a reasonable range, but Baccarani 
and Wordeman point out in |H] that neglecting thermal conductivity completely can lead to 
nonphysical results and mathematical instability. Although their work is directed to Si, their 
remarks should be equally valid for GaAs since the equations have a similar form in both cases. 
In our simulations, wc have chosen r = —2. 

2.3. The energy balance model 

The EBM is obtained as a simplification of the full HDM by neglecting the convective terms 
of the momentum balance equation (jS)). Additionally, the energy balance equation is 
simplified by the assumption that the time derivative of the mean electron energy duj/dt 
is small compared to the other terms and that the kinetic part in to can also be neglected, i.e. 



LO — —nkT , 
2 

or, if we take the two- valley structure of GaAs into account. 



(23) 




(24) 



The energy balance equation then becomes 



envE - V{nkTv) - V{-kVT) 



UJ — 



- ^nkTL 



(25) 
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and the momentum balance equation simplifies to the well-known current equation 

e -I e ^/nkT\ 

J = TpuE TpV . (26) 

m m \ e / 

The continuity equation and the Poisson equation remain of course unchanged in the EBM. 
Neglecting the time derivative of the current density is equivalent to the assumption that the 
electron momentum is able to adjust itself to a change in the electric field within a very short 
time. While this assumption is justified for relatively long-gated field effect transistors, it needs 
to be investigated for short-gate cases. 

Setting the electron temperature T equal to the (constant) temperature of the semiconductor 
material leads to the drift-diffusion model. Such a simplification is clearly not justified for 
the case studied in this paper (see also 0). 



3. NUMERICAL ASPECTS 
3.1. Discretization of the equations 

Today, many elaborate discretization methods are available for the DDM equations or EBM 
equations. The well-known Scharfetter-Gummcl method jJOl for the DDM makes use of the fact 
that the current density is a slowly varying quantity. The current equation is then solved exactly 
under the assumption of a constant current density over a discretization cell, which leads to 
an improved expression for the current density than it is given by simple central differences. It 
is therefore possible to implement physical arguments into the discretization method. Similar 
techniques have been worked out for the EBM ^T] . But due to the complexity of the HDM 
equations, no satisfactory discretization methods which include physical input are available 
for this case. 

Therefore, we developed a shock-capturing upwind discretization method, which has the 
advantage of being simple and reliable. For our purposes, it was sufficient to use a homogeneous 
mesh and a constant time step. But the method can be generalized to the non-homogeneous 
case. 

Time discretization is done for all equations by forward Euler differencing, i.e. the 
discretization scheme is fully explicit. The discretization is always written down only for the 
x-component of vectorial quantities in the sequel, since the corresponding expressions for y- 
components are then easy to derive. 

The constant timestep At used in our simulations was typically of the order of a few tenths 
of a femtosecond, and quantities at time T = tAt carry an upper integer index t. 

The rectangular simulated region of the MESFET gets discretized into x Ny rectangular 
cells Cij of equal size Ax x Ay — (Ix/N^) x (ly/Ny). Scalar quantities at timestep t like 
nl jjUjj jjTlj and where i — 1, ...N^ and j — 1, ...Ny are located at the center of the cells, 
whereas vectorial quantities like e.g. the electric field components E*^.i^i/2ji i j+1/2 '^^ 
velocity components ^^.^+1/2^1 ''^y i j+1/2 ^'"^ always calculated first at midpoints between the 
scalar quantities. 

If necessary, we can define intermediate values, e.g. E^-ij by 

Ex-i,j = —{Ex-i^l/2,j + Ex-i+l/2,j) , (27) 
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but a different definition applies to e.g. jx-.ij, as we shall see. 

The fundamental variables that we will have to compute at each timestep are riij, (pij, 
LOij (or Tij), jx-i+i/2,j and always respecting the imposed boundary conditions. All 

other variables used in the sequel should be considered as derived quantities. 

The momentum balance equation is discretized in the following way: 



k 



- -1"'i+l/2,j^x;i+l/2,j 

"(Px;i+l/2j'*^2:;i+l/2j ^ V^i-l l2,j'^x;i-\ I2,j) I (28) 

-(P^;i+l/2J^'!/;iJ+l/2 -pL;i+l/2j-l^!/;iJ-l/2)/Ay (29) 

~Pa;;i+l/2,j/'^p;i+l/2,j ' ("^0) 



where Pa;;i+i/2,j > and Py;ij+i/2 > and the same discretization strategy is applied to the 
y-component of the electron velocity. Prom the momentum density we obtain the new particle 
current density by 

■?x^+l/2J = -Pxti+l/2,j/"^i+l/2J ' (^1) 

and the momentum density at is extrapolated from neighbouring points in the direction 
of the electron flow x-component 

( 3 t+l _ 1 i+1 ^ . t+1 ^ n 

t+1 _ I 2^a:;i-l/2,i 2fx;i-Z/2,j> ■ ^x;i+l/2j ^ " (or,\ 
Px;i,i — 1 3^t+l l^t+1 \ . ^t+1 , n ) l^^^j 



and finally we get 



3j,«+^ - 1 • < 

2^x;i+l/2,j 2t^x;i+3/2,j) ' ^x;i+l/2,j ^ " 



<li:j=P^jKj/rr^!j > (33) 



V 



/2,j - jx;i+l/2,j/'^i+l/2,j/'^i+l/2,j ■ (34) 

We found that the purely heuristic choice 



"*+l/2J = \J<A+h3 (35) 

in the equations above improves the stability of our code. 

The electron temperature is related to the energy density by the relation uj\ - = ^nl jkT^j + 
\fn**jn\ j{v*^.^j +^^y?i,j) + /^Lii.j^-^rz, and can therefore be regarded as a dependent variable. 
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The upwind discretization is of the energy balance equation is given by 



At 









Oj': ■ 






1 

~Ax 


r* ■ • 

0a;;e,j+l/2j ~ 


jx;e,i-l/2,j) 


1 

Ax 


(jx;p,i+l/2,j ~ 


3x;p^i— 1/2, j) 


1 

Ax 


Ux;h,i+l/2,j ~ 


jx;h,i-l/2,j) 


1 


{fy;e,i,j+l/2 ~ 


jl;e,i,j-l/2) 


1 




•^y;Pi»J-i/2) 






1 


ijl;h,i,j+l/2 - 


fy;h,i,j-l/2) 



(36) 



where we have defined the energy currents 

n't — -,,*+! , ,t 



jx;e,i+l/2,j — '"x;i+l/2,j'^i+l/2,j ' (37) 
jl;p,i+l/2,j = ^ix^+l/2,j^/+l/2,j ' (^^8) 

and 

jx;h.i+l/2,j = ~'^i+l/2,ji'^i+l,j ^ "^ij)/ ■ 

Having obtained the new values for the mean electron energy, the transport parameters and 
energy-dependent masses are then also updated. 

The current continuity equation is discretized in a conservative way, using the particle 
current density j = nv 

- n* - 

~ ~ ~Ux;i+l/2,j ~ jx;i-l/2,j)/^^ 
-iJ'y■,^,J + l/2-Jr,^,J-l/2)/^y , (40) 

i.e. particles that leave cell in x-direction enter cell (i + and analogously for the 
y-direction; therefore, the total number of electrons inside the MESFET can only be changed 
at the boundary of the simulation region (mainly at the contacts). 

The Poisson equation, which is coupled to the hydrodynamic equations only through the 
particle density n, can be solved by any convenient method which is fast enough, since the 
computational effort should be kept as small as possible. Therefore, we used a multigrid method 
to perform this task. Fortunately, the Poisson equation has not to be solved at each timestep. 
Since the relaxation times of GaAs are of the order of some tenths of a picosecond, a timestep 
of about ten femtoseconds is fully sufficient for the update of the electric field. Fortunately, 
the stability of the discretization scheme is not affected that way, and allows an enormous 
reduction of the computational costs. 
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3.2. Boundary conditions 

For the basic quantities n, (j) and lo we imposed Dirichlet boundary conditions at each timestep 
at the contacts, e.g. the potential 4> at the source and the drain was fixed by the apphed voltage 

0Im = VsA- (41) 

Analogously, we assumed charge neutral contacts at the source and the drain, such that the 
charge carrier density was given there by the fixed doping density. The gate contact was 
modelled by assuming a Schottky barrier height of 0.8 V. For further details concerning the 
metal-semiconductor contact modelling we refer to standard textbooks . As far as the energy 
density is concerned, we imposed the boundary conditions directly on the electron temperature 
by assuming that the electron gas is in thermal equilibrium with the drain/source contacts. 

The artificial boundaries were modelled using von Neumann type boundary conditions. We 
present an explicit example in order to illustrate this point. We assume that the discrete values 
of the electron density at the MESFET boundary between the source and the gate (see Fig. 
5) are given by ...n^^i^ iy. The first index denotes the horizontal direction, whereas the 

second index starts with 1 at the top of the MESFET. Then, after each update of the density 
according to eq. (|^ . we enforce 

•■•"-(»2,i) = "-(^1,2), •■•"-(J2,2), (42) 

corresponding to the von Neumann condition that the normal component of the electron 
density vanishes at the specified boundaries. 

As mentioned above, we used a multigrid algorithm to calculate the electric potential. Also 
there, the mixed Dirichlet/ von Neumann boundary conditions were imposed on the subgrids 
at each intermediate step of the calculations. A FORTRAN90 program which calculates the 
potential by a multigrid algorithm can be obtained from the first author's address. 

4. SIMULATION RESULTS 

The GaAs MESFET structure used in our simulation is shown in Fig. 5. The structure consists 
of a 0.1 /xm-thick active layer with a doping concentration of Nd = 2 • lO^^cm""^ on a 0.3 /im 
buffer layer [Nd = lO^^'cm"^). The doping profile is abrupt between the two layers, the lattice 
temperature is = 300 K. For steady-state results, we used long simulation times of 30 ps 
such that the steady state was de facto reached. The length of the drain and source contacts 
is 0.5 /im, the gate-source separation 0.5 /xm, the gate-drain separation is 1.0 /xm and the gate 
length is 0.8 /^m. The Schottky barrier height is assumed to be 0.8 V. 

In order to obtain stable and physically meaningful results, values like Ax — Ay — 6.1 
nm for a grid size of 537 x 65 were typical values used in the simulations. The rather large 
mesh was manageable due to the effective multigrid algorithm used for the solution of the 
Poisson equation. We found that the mesh size used in |F| was too coarse for accurate results, 
although the authors improve accuracy and convergence speed of their calculations by using a 
non-homogeneous grid. In fact, a non-homogeneous grid necessitates additional calculational 
costs which reduce the speed of the simulation, and the timesteps also depend on the size of 
the smallest cells. Furthermore, the Poisson equation was solved by a conventional successive 
overrelaxation method. 
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Figure 5. The MESFET geometry. 
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Figure 6. Electron temperature inside the MESFET. 



We present results for a gate-source bias Vgs = V and a drain-source bias Vds = 3 V. 

Fig. 6 shows the electron temperature inside the MESFET for the stationary case. 

Fig. 7 and 8 show the electron velocity and the electron temperature along the channel of 
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Figure 7. Electron velocity inside the MESFET for a fine grid (solid line) and a coarse grid (dash- 
dotted). 



the MESFET (0.077 below the contacts). Due to the high temperature which is reached 
under the gate, the electron mobiUty is strongly lowered in this region. The electrons, which 
overshoot under the gate, are therefore deaccelerated abruptly to a lower velocity in the high- 
temperature region. The results obtained from an energy balance calculation are in good 
quantitative agreement, the differences are mostly pronounced in the region where the electron 
velocity is high, as it is expected from the different treatment of the energy density in the HDM 
and EBM. 

In the very close vicinity of the gate, it was necessary to reduce artificially the electric field 
or the velocity of the electrons in order to stabilize our code. We checked that this does not 
strongly affect the simulations results outside this region due to the very low density of the 
electron gas near the contact; a similar procedure was also necessary in 

The velocity and temperature curves are in fact very similar to those of one-dimensional 
simulations of ballistic diodes (n+ — n — n^-structures), which were used as simplified models 
for FET channels. 

It is interesting to observe in Fig. 7 that our simulation results for a coarse grid (grid size 
217 X 33) are very close to those presented in Fig. 5 in [S], where a non-uniform mesh of 
typical size 141 x 35 was used for a similar MESFET geometry. This may be considered as a 
confirmation of the results given in for a relatively coars grid, but demonstrates the fact that 
stability of the code does not automatically imply accuracy of the results, and an investigation 
of the dependence of the results on the grid resolution is indispensable. The channel below the 
gate (see also Fig. 10), which is the most interesting region in the MESFET, is relatively small 
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Figure 8. Electron temperature inside the MESFET along the channel. The dash-dotted line shows 

the EBM result. 

and must be resolved sufficiently. 

Due to the strong heating of the electron gas in the channel region, most electrons are excited 
to the L-band. This leads to a cooling effect of the electron gas, since the excitation energy 
is missing in the thermal energy balance. The dash-dotted curve in Fig. 9 shows the electron 
temperature that would be obtained from the HDM if the energy term /?L(ZU)Ai?rL in eq. 
were neglected. The energy AEtl = 0.29 eV which is necessary to excite an electron to the 
upper conduction band corresponds to a temperature difference AT = 2Ai?rL/3A: ^ 2000if ; 
the observed error is of the same size. 

Finally, Fig. 10 shows a surface plot of the electron density inside the device. Clearly visible 
is the MESFET channel under the gate. 



The feasibility of two-dimensional hydrodynamic simulations is demonstrated for the case of 
a GaAs MESFET structure. Although the single-gas hydrodynamic model is superior to the 
drift-diffusion or energy balance model, it is desirable to direct the efforts of future research 
in the direction of multi-valley hydrodynamic models. Models like the EBM will no longer be 
adequate for the physical description of high-speed submicron devices in the near future. 

It is obvious that future attempts to model submicron devices will face many more problems 
which have not been touched in this paper. One difficulty is the fact that the components 
of semiconductor devices are often of very different size and material composition. This 
necessitates the use of adaptive discretization grids or the hybridization of different numerical 
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Figure 9. Electron temperature inside the MESFET along the channel for the 'correct' model (solid 
line) and the 'wrong' model where the excitation of electrons to the upper conduction band is not 

taken into account properly. 
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methods, but both strategics arc hampered by severe problems hkc mimcrical instabihties 
or huge computational efforts for realistic simulations. Another class of problems is due to 
the fact that the physics of semiconductor materials is very complex, and therefore hard to 
implement such that the physical behavior of the device is satisfactorily described. Attempts 
to describe quantum effects in device modelling should be considered as tentative in the case 
of heterostructure devices. There is no optimal solution for these problems, and the numerical 
and physical models have to be adapted to the problem under study. We hope that our detailed 
description of a hydrodynamic simulation may serve also as a help for researchers entering the 
field of hydrodynamic semiconductor device modelling. 
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